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We obtain an analytical expression for the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction 
J in electron or hole doped graphene for linear Dirac bands. The results agree very well with the 
numerical calculations for the full tight-binding band structure in the regime where the linear band 
structure is valid. The analytical result, expressed in terms of the Meijer G-function, consists of a 
product of two oscillatory terms, one coming from the interference between the two Dirac cones and 
the second coming from the finite size of the Fermi surface. For large distances, the Meijer G-function 
behaves as a sinusoidal term, leading to the result J ~ R~ 2 kF sin(2fcF-R){l + cos [{K — K').R]} for 
moments located on the same sublattice. The R~ 2 dependence, which is the same for the standard 
two-dimensional electron gas, is universal irrespective of the sublattice location and the distance 
direction of the two moments except when fci? = (undoped case), where it reverts to the R~ 3 
dependence. These results correct several inconsistencies found in the literature. 

PACS numbers: 75.30.Hx; 75.10.Lp; 75.20.Hr 



The Ruderman-Kittel-Kasuya-Yosida (RKKY) 
interaction,— which measures the coupling between 
two magnetic moments mediated by a background of 
electrons, is an important characteristic of the electron 
system. It has been extensively studied for the electron 
gas in one, two, or three dimensions. Even though 
graphene is a two-dimensional (2D) system, there are 
two important differences from the standard 2D electron 
gas, viz., the linear band structure and the existence of 
two Dirac cones in the Brillouin zone. This leads to the 
two characteristic momenta, the Fermi momentum fcp 
and the momentum difference of the two Dirac points 
K — K', both of which produce oscillatory factors of 
their own, leading to unusual features not found in the 
standard 2D electron gas, e.g, the beating of the RKKY 
interaction that can be controlled by a gate voltage. 

Although the RKKY interaction in graphene has al- 
ready been studied by many authors^— the results differ 
from one another, even in the long-distance behavior, 
and the interference term from the Dirac cones is often 
missing in the results. In this paper, we derive analyti- 
cal expressions for the RKKY interaction for the linear 
band model, extending our earlier work for the undoped 
case.— The results are compared to the same for the tight- 
binding model which we also calculate from the numerical 
evaluation of the lattice Green's function. The linear- 
band and the tight-binding results agree quite well in the 
cases where fc^ lies in the linear regime (\Ep\ < t/2>). We 
find that the analytical results in the linear band approx- 
imation may be expressed as a product of the J for the 
undoped case, which is obviously independent of fep, and 
a new factor that depends on Uf and goes to one in the 
limit of kp — > 0, so that the results for the undoped case 
are correctly reproduced. The analytical results, sum- 
marized in Table I, are expressed in terms of the Meijer 
G-function, whose long distance behavior is sinusoidal. 

Model and the method - We consider the nearest- 
neighbor tight-binding Hamiltonian for the 7r-electrons 
in graphene including the contact interaction with two 
localized magnetic moments 



where Ho = —t^2uj) a c i<T C j°- + H- c - ^ s tne tight-binding 
Hamiltonian, (ij) denotes summation over distinct pairs 
of nearest neighbors, t ~ 2.56 eV,— a is the spin index, 
and the interaction term between the localized spins S p 
and the itinerant electron spins s p is given by "Hint = 
— A(5i ■ si + S2 ■ 82). In the linear response theory, the 
interaction energy may be written in the Heisenberg form 



E(R) — J a p(R)Si ■ D2, 



(2) 



where the sublattice indices and the positions 
of the two moments are (a, 0) and (j3,R), 
J aj3 {R) = 4- 1 X 2 h 2 x a i3(0,R) is the RKKY interaction 
and Xaf){ r , r ') = Sn a {r)/SV/3(r') is the susceptibility. 
Note that R denotes the position of the atom and not 
the position of the cell in which it is located; they differ 
by the basis vector of the atom in the unit cell. 

The susceptibility can in turn be computed from the 
unperturbed Green's function^ - — 



dE Im[G° a J0,R,E)G°s a (R,0,E)] 



. (3) 

Below we evaluate this integral analytically for the linear 
bands and numerically for the tight-binding bands by 
direct integration 

G%(R,0,E) = ^~J d 2 ke ik - R G° a0 (k,E), (4) 



of the momentum-space Green's function 

E + in + Uk 



G%{k,E) 



(E + ir,y-\f(k)f 



(5) 
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HeieHk = I f*ffc) J i s tne § ra P nene tight-binding 

Hamiltonian in the momentum space and the Bloch sum 
f(k) = -t (e* dl + e ik - d2 + e lk - d3 ), where di, d 2 and d 3 



1~L — 7"^o ~t~ 



(i) 
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TABLE I: Summary of the RKKY interaction in graphene for both the doped (fop 7^ 0) and the undoped case (kp = 0), 
given as a product of the terms: J a p = otcotDOtF- The long-distance behavior is obtained by replacing of with a' F . Here 
C = 9X 2 h 2 /(256nt), C = -\' 2 V 2 m* /8(Nnh) 2 , x D = (K - K') ■ R, x F = k F R, 9r is the angle of the position vector R made 
with K' — K direction, where K' and K are two adjacent Dirac points in the Brillouin zone. The results for the standard 
two-dimensional electron gas (2DEG) 1( ^ are also shown, where we have rederived the long-distance behavior. 
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are the three nearest-neighbor position vectors. Note 
from expressions following Eq. © that the Friedel 
oscillations^ Sn a (r) in the charge density induced by a 
(5-function potential is proportional to J a p as well. 

Moments on the same sublattice - Using methods dis- 
cussed in our previous work^, the Green's functions as 
well as the susceptibility can be evaluated both for the 
linear-band approximation and for the full tight-binding 
bands. For the linear-band case and for moments on the 
same sublattice, the result is 

Xaa{0, R) = Iaa(R) x {1 + cos[(if - K') ■ R}}, (6) 

where 



Iaa(R) = -- [ F dElm [g AA (R,E)f 

™ J-oo 



(7) 



and using the integral tables^ along with fj, = 2, v = 0, 
and the new variable x = z 2 (fci?i?) -2 , the result is 

dz z 2 J (z)Y (z) = FU / dx x-^ 2 x 

° 3 I ( 10 ) 

G M (l,f,l k2pR2x 



k F R 
2^ 



M(k F R), 



( 13 

where M(k F R) = , 2 \ 2 -1 

hand notation for the Meijer G-function. Using Eqs. (JB]), 



kpR 2 I is a short- 



(|8]), and (|T0|) . we arrive at our desired result, valid for the 
moments on the same sublattice and for the linear bands, 
viz., 

J AA {R) = J° AA (R) [1 + ^£^ M (k F R)], (11) 



g AA {R,E) = -2TrEVp 2 nszK (-iER/v F ), K is the 
modified Bessel function of the second kind, v F — 3ta/2 
is the Fermi velocity, a is the bond length and Qbz is the 
area of the Brillouin Zone. Now we split the integral in 
Eq. ([7]) into two parts, viz., = f®^ + f® F , where 
the first term accounts for the valance electrons (undoped 
case) and the second for the conduction electrons, so that 



Iaa(R) = ^—R-^Io- 



k F R 



dz z 2 Jo(z)Y (z)}, (8) 



where I = - J °° dy y 2 J (y)Y (y) = -1/16, 5 y = 
—ER/v F for the valance band (E < 0), z = ER/vf 
for the conduction band (E > 0) and Jo and Y Q are the 
Bessel and Neumann functions with real arguments and 
kp is the Fermi momentum. 

The remaining integral in Eq. (j5]) may be expressed 
in terms of the Meijer G-functions. The product of the 
Bessel and the Neumann functions can be written as 



z»J v {z)Y v {z) 



r-V2 



"1,3 



2 ' 2 



2 



(9) 



where J\ A {R) = -C x {R/a)~ 3 {l + cos[{K - K') ■ 
R]} is the undoped exchange interaction with C = 
9A 2 fi, 2 /(2567rf). The only approximation used here was 
to extend the linearity of the Dirac bands to infinity (in- 
finite momentum cutoff); however, this approximation is 
in good agreement with the numerical full-band tight- 
binding calculations, both for the undoped case^ and for 
the doped case if kp is small [Fig. JT}]. 

Note that in the expression for the RKKY interaction 
Eq. (jlip the Fermi momentum term in the square bracket 
depends only on the magnitude of the distance R, while 
the Dirac cone term J\ A (R) depends on its direction as 
well, which makes the interaction direction dependent. 
Here K and K' are any two adjacent Dirac points in the 
Brillouin zone. It is easy to see that while the oscillatory 
factor 1 + cos((K — K') ■ R) repeats in triplets as 2, 1/2, 
1/2, ... with distance R along the zigzag direction, it 
is always two for the armchair direction, so that J AA 
changes smoothly along the armchair direction but not 
for the zigzag direction [Fig. @]. 

One is often interested in the long-distance behavior 
of the RKKY interaction and this may be obtained from 
the asymptotic behavior of the Meijer G-function M(x). 
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FIG. 1: (Color online) RKKY interaction Jaa obtained for 
the tight-binding bands (solid line), compared to the RKKY 
expression involving the Meijer G-function Eq. (red dots) 
and its long-distance limit Eq. (|14p (dashed line). 



We find using standard tables^ that 

4x 2 [l-37-31n(a;/2)] 

— 2cos(2x) — 8a;sin(2a;) 



lim M(x) 



lim M(x) 



(12) 



(13) 



where 7 ss 0.577 is the Euler-Mascheroni constant. These 
functions are plotted in Fig. ([3]). We note that with this 
asymptotic dependence, the square bracket in Eq. (fTTj) 
becomes one for fc F = 0, so that the RKKY interaction 
becomes the same as for the undoped case as it must. 

Contrary to J\ A (R) that always shows ferromagnetic 
coupling for the moments on the same sublattices owing 
to the particle-hole symmetry^, the oscillatory behavior 
of M(kpR) leads to the oscillations of Jaa between ferro- 
magnetic and anti-ferromagnetic interactions. From Eqs. 
(fT3| and (fTT|) . we obtain the long-distance behavior 

lim Jaa(R) = J^ j4 (i?)7r" 1 [2cos(2a;F)+8a; F sin(2a; F )], 

A' y R — > OO 

(14) 

where xp — kpR. Note that the distance dependence is 
R~ 2 if kp is non-zero, i. e., the same as for the ordinary 
2D electron gas^i. If kp — 0, the RKKY interaction 
reverts to the undoped case as seen from Eqs. (fTTj) and 
([T2"T) . so that the distance dependence is now R~ 3 . 

It is worth mentioning that the correct result for 
kpR >• 1 can only be found by evaluating the 
Meijer G-function for large arguments and not just 
by replacing the Bcssel functions in Eq. (|10|) by 
their large-argument (z 3> \v 2 — 1/4| ) limits, viz., 
J v (z) « 2 1 /2( 7rz )-V2 cos ( z _ ly7T /2 - tt/4) and Y v {z) w 
2 1/2 (7rz) _1/2 sin(z - ^7r/2 - 7r/4) . The latter approach 
happens to lead to the same functional form as in Eq. 
([T^]) but with incorrect coefficients because of the error 
made in the small kpR contribution to the integral in Eq. 
(fTUfl . However, there is no such problem if kpR <C 1 and 
the short-range results can be found either way to yield 

32(k F R) 3 



J AA {R) = J° aa (R){1 + - 



9tt 



-[l-3 7 - 



•31n(fe F i?/2)]}. 

(15) 
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FIG. 2: (Color online) RKKY interaction for several cases. 
Black solid lines are the numerical results for the full tight- 
binding band structure, while the red lines indicate the ana- 
lytical results, Eqs. and (I19[) . for the linear bands. 




FIG. 3: (Color online) The Meijer G-function as a function of 
k F R (black solid line) versus its asymptotic expansions (Eqs. 
Upland 113)1, 



Moments on different sublattices - For moments lo- 
cated on two different sublattices, we proceed as before 
to obtain the susceptibility^ 

X ab (0, R) = I AB (R) x { 1 + cos[(K - K') ■ R + tt - 26 R ] }, 

(16) 

where I AB {R) = f J-L dE Im [9ab{R,E)} 2 , and 
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FIG. 4: (Color online) Switching the exchange interaction 
between ferro and antiferro by changing the carrier density in 
graphene with gate voltage. The larger the distance between 
the moments, the earlier is the switching, which is controlled 
by k F R- 
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FIG. 5: (Color online) Beating pattern of the RKKY inter- 
action for bond-centered moments separated along the zigzag 
direction. 



cjab{R,E) = -2itEv F 2 Q,' B 1 z K 1 {-iER/v F ), Expanding 
the modified Bessel function K\, the integral becomes 



Iab(R) 



3e>-3 



lir^R 



n 2 BZ v F 



k F R 



[I + / dz z^J 1 {z)Y 1 {z)], (17) 



where Iq = — / °° dy y 2 Ji(y)Yi(y) = 3/lG 2 -^ is the contri- 
bution from the undoped part and the remaining integral 
can again be expressed in terms of the Meijer G-function 



k F R. 



dz z 2 J 1 (z)Y 1 (z) = -'^M'{k F R). (18) 



/o 

This leads to the final result 



Jab {R) = J AB (R) [1 - ^M'(k F R% (19) 



where M'(k F R) = G 
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2> 2 



1,2,0,^ 



k 2 F R 2 



undoped exchange interaction is J AB (R) 
(R/a)- 3 {l + cos[{K -K')-R + n- 26 R ]}. 



and the 

= 3C x 



It has been demonstrated that the dopant carrier con- 
centration in graphene can be controlled by a gate voltage 
or chemical doping. 12 This raises the interesting possibil- 
ity of switching the magnetic interaction between ferro 
and antiferro. This is illustrated in Fig. [4j where the 
exchange interactions were evaluated using the full tight- 
binding bands as a function of k F and the carrier density 
is given by n — k F /ir in the linear-band region. 

Bond moments and the beating of J - For moments 
located on the bond center, the interaction is of the form 
Hint = —\S -J2 P $p, where the summation is over the two 
adjacent atoms. The exchange interaction becomes the 
sum of the site interactions: Jbond = 2 Jaa+Jba+Jab = 
Sn-^-CiR/a)- 3 x {2cos(2fc F i?) - 3 cos(2k F R) cos[(K - 
K')-R]~Ak F R sm(2k F R)cos[(K-K')-R}}. The resulting 
beating pattern of the RKKY interaction is shown in Fig. 
(O, which can be controlled by the gate voltage. 

Finally, we note that Table I is valid both for elec- 
trons and holes because of the particle-hole symmetry, 2 
Mathematically, this follows from the fact that the net 
contribution to the susceptibility from a symmetric range 
of energy is zero as may be seen by taking the integral 
in Eq. ([3]) from — e to e and by using the symmetry: 
G° afj (0, R, E) = G° pa (R, 0, E) and the fact that the prod- 
uct Im G° (3 (0, R, E) x Re G° af3 (0, R, E) is an odd function 
of energyiii^ 

In summary, we provided analytical results for the 
RKKY interaction in graphene in the linear-band ap- 
proximation and showed that these results agree with the 
numerical results obtained for the tight-binding bands if 
the Fermi momentum is small. The presence of the two 
characteristic momenta, viz., the Dirac cone momentum 
K — K' and the Fermi momentum k F , leads to the un- 
usual oscillatory features in graphene. 
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